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The critical behavior of the 0(2) model on dilute Levy graphs built on a 2D square lattice is 
analyzed. Different qualitative cases are probed, varying the exponent p governing the dependence 
on the distance of the connectivity probability distribution. The mean-field regime, as well as the 
long-range and short-range non-mean-field regimes are investigated by means of high-performance 
parallel Monte-Carlo numerical simulations running on CPU's. The relationship between the long- 
range p exponent and the effective dimension of an equivalent short-range system with the same 
critical behavior is investigated. Evidence is provided for the effective short-range dimension to 
coincide with the spectral dimension of the Levy graph for the XY model. 



I. INTRODUCTION: INTERACTING MODELS 
IN COMPLEX NETWORKS 



The study of interacting systems defined in complex, 
non-regular structures is interesting from at least three 
points of view. Statistical mechanical models in graphs 
are used for the description of phenomena in differ- 
ent fields, among which one can cite: stock market 
dynamics correlations in bird floc king ,^ avalanches in 
brain activity® or biological networks!^ Furthermore, 
there is a theoretical interest per se in the study of crit- 
icality in complex networks. The theory of critical phe- 
nomena establishes that the critical properties of systems 
interacting in a d-dimensional lattice only depend on the 
symmetries of the interaction and on the dimensionality 
d. On the other hand, when the topology of the inter- 
action is more complicated, e.g., translational invariancc 
is lost and symmetries of the lattice are broken, the de- 
pendence of criticality on the structural properties of the 
interacting graph is not known in general, although this 
topic has been subject of interest since more than two 
decades and several results are available for particular 
models P 

A particularly clear case is the spherical model, which 
has been proved to be equivalent to the n — > oo limit 
of the 0(n) model when both models are defined on a 
latticed On general graphs the full equivalence does not 
hold anymore, though the critical behavior of the spher- 
ical and 0(n — > oo) still do coincide.^ The critical prop- 
erties of the spherical model on a general graph are ex- 
actly knowrP^ anc j are such that the universality class 
of the transition only depends on a single quantity: the 
spectral dimension of the graph, d, defined in terms of 
the low-frequency spectral density of its adjacency ma- 
trix q{lo) ~ w d / 2_1 . This quantity is also related with 
the probability of self-return of a random walker in the 
graph, and determines the infrared div ergenc es of a Gaus- 
sian field theory defined on the graphP2LG3 Remarkably, 
the functional dependence of the spherical model critical 
quantities on a graph with spectral dimension d turns 
out to be the same as the ones of a short-range model in 
a hyper-cubic lattice with Euclidean dimension d. This 



analogy provides a suggestive, physical sense to the non- 
integer dimensions appearing in the context of the theory 
of critical phenomena. 

The spectral dimension also plays a role in the XY 
model criticality, which was proved^ to exhibit sponta- 
neous magnetization in the ordered phase in a graph of 
spectral dimension d > 2, and absence^ of spontaneous 
magnetization for d < 2. The latter phenomenon being 
well known in the two-dimensional (2D) XY modelp^HII] 
which is a particular case of this result. 

Further numerical and analytical results for the criti- 
cality of other particular models in graphs are available 
(for a review see Ref. [5]). The Ising model was first 
studied in Small- World networks p^HiS] jnBarabasi- Albert 
networks p2H2U anc j on general graphs J2S12ZI w here it was 
found that the universality class depends on the diver- 
gence or finiteness of the second and fourth moments of 
the degree distribution. In this way, three different crit- 
ical regimes may be discriminated: (i) absence of phase 
transition, when both second and fourth moments di- 
verge; (ii) a non-mean-field second order transition, when 
the second moment is finite; and (iii) a mean-field second- 
order transition, when both moments are finite. 

Studies of the Ising model in scale-free netwo rks^ 
and in correlated growing-random network d 29 * 30 ! were 
also performed. In the latter case a phase transi- 
tion was found of the Kosterlitz-Thouless (KT) uni- 
versality class, different from the mean-field nature of 
the transition found in (uncorrelated) scale-free net- 
works. This difference was argued to have its origin in 
the sign of the degree-degree correlations (ass ort at i vity- 
disassortativity) of both types of networks ! 29 * 31 ! The 
Potts model has also been investigated, 32 34 finding an 
infinite-order transition for a divergent second moment 
of the degree distribution. 

Eventually, also the 0(2) XY model has been analyzed. 
In the ID Small- world network, it was arguecP^to exhibit 
long-range order for arbitrarily low values of the rewiring 
probability (like in the Ising case). For uncorrelatecP and 
correlated^ scale-free networks an order-disorder transi- 
tion is observed for a sufficiently large value of the degree 
distribution exponent. Interestingly, as it happens in the 
Ising case, in the correlated scale-free network the transi- 
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tion is non-mean-field, unlike the uncorrelated case. This 
difference is again ascribed to the different nature of the 
degree-degree correlations in both kinds of graphs. 4 

Another piece of the puzzle is provided by the numeri- 
cal work carried out by Yang et al., 36 in which the critical 
behavior of the XY model is studied in uncorrelated and 
correlated random (rather than scale-free) graphs. In the 
first case, i.e., the Erdos-Renyi graph, the transition is 
found to be of the mean-field type, while in a randomly 
growing network it is claimed that the occurring transi- 
tion belongs to the KT universality class. 

Despite numerous results in this field, a unified picture 
of critical phenomena in graphs is still lacking. For in- 
stance, it is not clear under what conditions a relation can 
be established between criticality in graphs with spectral 
dimension d and short-range models in d dimensional lat- 
tices, nor what is the relation with the conjectured influ- 
ence of dissasortativity on criticality. 

Considering the field of disordered systems as well, 
spin-glass models defined in complex graphs have been 
recently investigated, on the so-called Levy or dilute 
lattice^ that we will use in the present work. It is a 
graph whose sites are connected with a probability de- 
caying as a power p of the distance between them, as 
computed in a given, auxiliary, nearest-neighbor lattice. 
While for large enough p one recovers the auxiliary lat- 
tice, the p = Levy graph limit corresponds to the ran- 
dom Erdos-Renyi graph, such that the zN/2 bonds {z 
being the coordination number of the lattice) are cho- 
sen at random from the set of all N(N — l)/2 possible 
bonds. The Ising spin-glass model in the Levy graph 
was introduced in Ref. EZI as a dilute variant of the fully 
connected Ising spin glass with long-range interactions in 
one dimensionpSHU to allow for more efficient, i.e., 0[N] 
computation of observables in a graph with N sites. Like 
the fully connected long-range spin glass, the dilute one 
is such that, varying the power p, one actually acts as 
if varying the dimension of a L>-dimensional short-range 
lattice model, equivalent -from th e criti cal behavior point 
of view- to the long-range modelf^EI] 

Starting from the field-theoretic representation in the 
free theory limit, valid both for ordered and disordered 
systems, an equivalence between p and D can be conjec- 
tured to btP 
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p £ (d : 2 + d] 



D = 

p — d 

D=d p>2 + d 



(1) 



where d is the dimension of the auxiliary lattice of the 
dilute model. This can be improved as^ 



D 



(2 - Vsr)d 



P- 



(2) 



where f? sr is the anomalous scaling exponent of the space 
correlation function in the Z?-dimcnsional short-range 
counterpart. The above relationship does not depend on 
the specific symmetries of the system, nor on the pres- 
ence of any long-range order at all (as in the quenched 
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TABLE I: Summary of the different dimensions considered. 



disordered case) . What can change is the range of values 
of p determining the universality class at which the model 
belongs. In particular, for fully connected systems with 
(ordered or disordered) long range interactions decaying 
with the /5-th power of the distance in a d-dimensional 
hyper-cubic lattice, three regimes can be identified: 

• d < p < p m {{d), in which the system undergoes a 
mean-field transition; 

• Pmi(d) < p < 2 + d, in which infra-red divergences 
take place, to be dealt with renormalization group; 

• p > 2 + d where the critical behavior is short-range- 
like. 

The value of p m f(d) depends on the specific theory and 
its symmetries, thus being different in ordered^ (p m f = 
3d/2) and disordered^ (p mf = 4d/3) systems. 

Critical exponents are functions of p as well, as, e.g., 
7] p = 2 — p + d for any p (the r\ long-range expo- 
nent is not renormalized) and v p = l/(p — d), in the 
mean-field regime. Also these expression s are formally 
the same both in ordered^ and disordere d 39 * 42 ! systems, 
whereas different is their dominion in p and the renor- 
malized expression for p > p m [(d). The prediction for 
the rjp exponent has been compared with the outcome 
of numerical simulations in the case of the long-range 
Ising ferromagnetic On the other hand, Eq. pi has 
been carefully tested in ID Levy Ising spin-glasses for 
P > pmf(l) = 4/3, verifying that the equivalent short- 
range critical behaviors are actually consistent both for 
D = 3 (p = 1.792) and for D = 4 (p = 1.58)P 

To make reading more fluid and avoid notation ambi- 
guities, in Tab. [I] we summarize the various dimensions 
we refer to in this work. 

The possibility yielded by Levy lattices of changing the 
effective dimensionality, freely choosing the universality 
class of the model without compromising the complexity, 
is useful to approach different problems: the applicabil- 
ity of the replica symmetry breaki ng th eory in and out 
of the spin-glass mean-field regime EOill the existence of 
the Almeida-Thouless critical li ne ab ove the spin-glass 
upper critical dimension in Ishu j 48 ! 49 -! and Heisenberg^ 
systems, the criticality of the 3-spin spin glassp^ and the 
low temperature behavior in Heisenberg spin-glasses (in- 
cluding the spin-chirality decouplings 1 ^) as well as in 
0(m — > oo) spin-glasses PS 

There are some aspects that should be clarified also in 
this context. An elementary question is about the equiv- 
alence of the spectral dimension d of the graph and the 
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short-range equivalent dimension D, cf. Eq. (TJ. A rig- 
orous derivation of the critical properties of long-range 
dilute models as a function of the power p is still lack- 
ing, and also an argument stating under what conditions 
they are equivalent to the fully connected case with the 
same value of the power p. In what follows we will clarify 
some of the mentioned issues in the 0(2) case defined in a 
dilute 2D graph with open and periodic boundary condi- 
tions. Along with theoretical arguments, we shall present 
the outcome of numerical simulations run on Graphics 
Processing Units (GPUs) with a ad hoc optimized code, 
whose dynamics is based on the Metropolis, Parallel Tem- 
pering and Over-relaxation algorithms, suited to study 
continuous spin models interacting in graphs with arbi- 
trary topology. 

The paper is organized as follows: in Sec. [ll]we pro- 
vide a theoretical argument to support the existence of 
three intervals of the power p, relative to three different 
universality classes of the 0(2) model defined on a dilutc- 
2D graph with power p and we yield numerical evidence 
in support of the fact that D(p) is, actually, the spectral 
dimension of the graph with power p. In Sec. |III| nu- 
merical methods are explained and the implementation 
of our GPU-based algorithm and its performance are re- 
ported. We present numerical results in Sec. 
conclusions in Sec. |V] 
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FIG. 1: Probability distribution of the degree of connectivity 
of dilute 2D graphs with N = 256 2 = 2 14 nodes for three val- 
ues of the decay exponent p — 1/2, 10/3 and 14/3 of the link 
probability, one for each critical regime. Results are reported 
for both periodic (open symbols) and free (closed symbols) 
boundary conditions, p — 14/3 is in the short-range regime, 
p = 10/3 in the non-mean-field regime and p = 1/2 in the 
mean-field regime. In the latter case the Poisson distribution 
with average 4 is displayed for comparison. 
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II. CRITIC ALITY OF THE XY MODEL IN THE 
2D LEVY GRAPH 

We are concerned with the ferromagnetic 0(n) model, 
defined by the Hamiltonian: 

JV 

H = ~ / ! Jij s * • s j> ( 3 ) 

i<j=i 

where denotes the dynamic variable on the i th site 
of the graph, an n-dimensional vector with unit modu- 
lus, and the product is a n-dimensional Euclidean scalar 
product, the XY model being the n = 2 case. The values 
of the adjacency matrix J,^ of the graph can be either 
(no connection) or 1 (the interaction strength). We will 
study the dilute (Levy) graph, for which two sites i and 
j are connected (i.e., the element of the Jy matrix is 1) 
with a probability 

VpVii) = §|ri-ri|-' (4) 
Z = 

r 

and such that the total number of bonds is independent 
from p and equal to 2N. In Eq. Q, the vector cor- 
responds to the position of site i on a square lattice and 
the probability is normalized summing over the set of 
all possible distances between the sites of the 2D lat- 
tice. Operatively, the set of possible distances on lattices 



of linear size L depends on the boundary conditions cho- 
sen for the numerical simulation being periodic (PBC) or 
free (FBC). The maximum distance r max will be [L/2]y/2 
for PBC or Ly/2 for FBC. In the present paper the re- 
sults shown are for PBC, but the outcome of the nu- 
merical analysis presented in Sec |IV| remains unchanged 
when FBC are implemented, despite the differences with 
graphs with PBC in the p > case. In Fig. [I] we show 
the degree distribution of the dilute 2D graph (with free 
and PBC) for different values of p. While the p — > oo 2D 
limit of the Levy graph (i.e., the square lattice) exhibits 
a delta function S(d — 4), the p — limit corresponds 
to the Erdos-Renyi graph with degree distribution given 
by a Poisson distribution with average degree equal to 4. 
The latter case is independent from the kind of boundary 
conditions. 



1. A dimensional argument 

We now discuss the criticality of the model. Let us 
first consider the d-dimensional fully connected version 
of our model where each site is connected with any other 
site and the interaction strength decays with a power-law 
J(r) ~ \r\~ p of the distance in a c£-dimensional lattice. 
Following Ref. |44j we consider the following effective 
Ginzburg-Landau Hamiltonian for the long-range model, 
a scalar (j) 4 Sine-Gordon theory, 

K = L d J ^(^+m 2 )|0(q)| 2 + A J d<W 4 (r) (5) 
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where L is the linear size of the system, q is the momen- 
tum space index, m is the mass of the theory, (f> is the 
Fourier Transform of the scalar field, and A is the cou- 
pling strength. The long-range exponent ip is such that 
the Fourier Transform of the interaction J(r) goes like 



D 



J{q)=L- d ' 2 J d d v.J{v)e v 



•r-V> 



for low 
(J(r) ~ 



|q|. In the long-range fully connected case 
), it holds ip = p—d, and there is a divergence 
for p = d, the point at which the number of links, pro- 
portional to J(0), diverges in the thermodynamic limit. 

A dimensional analysis of the Hamiltonian, Eq. ([5]), 
see, e.g., Refs. 37 39 42 44, shows that the dimension 
of the coupling constant A is larger than zero whenever 
p > Zd/2. Below this point, A is an irrelevant variable at 
criticality: the system critical behavior is correctly de- 
scribed by a (mean-field) free theory. For p > 2 + d, 
the short-range lattice contribution to the propagator, 
q 2 , will take over the long-range contribution and the 
Ginzburg-Landau Hamiltonian will correspond to the one 
of a d-dimensional short-range model. One cold argue 
that these arguments still remain valid in the case of the 
dilute Levy lattice. We will motivate better in the follow- 
ing section such analogy (see also subsection II 2 1 . Even- 
tually, in the fully connected model, an ultra-extensive 
regime occurs for p < d, with diverging energy. This 
is not present in the dilute model, since the number of 
bonds of the graph is constant. Collecting all the above 
considerations, we summarize the following dependence 
of the criticality of the dilute XY model on the exponent 
p (see Fig. §: 

1. p > p sl - = 2 + d: the model behaves as its 
short-range version in d dimensions, for what con- 
cerns criticality, thus belonging to the Kosterlitz- 
Thouless (KT) universality for d = 2. This regime 
will be called the short-range regime in the follow- 
ing. 

2. p e (p m f : p sr ), with p m f = 3d/2: the system will 
present a transition different from a KT transition, 
with exponents different from the mean field-ones, 
this regime will be denoted as the (long-range ) non- 
mean field regime. 

3- p < pmf > the system belongs to the mean-field uni- 
versality class, i.e., its critical properties are the 
corresponding to a free Gaussian theory in dimen- 
sion 4. We will denote this regime as the mean-field 
regime. 



2. Spectral dimension 

Using Eq. 0, and naming D u the upper critical di- 
mension (D u = 4 for the present model) one checks that 
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FIG. 2: Three domains of p, relative to the three universality 
classes of the dilute XY model on Levy graph. The arrows 
point out at the values of p at which simulations have been 
carried out: 0, 14/6, 17/6, 20/6, 28/6. 



the regimes above introduced are in correspondence with 
the three regimes of an equivalent D-dimensional lat- 
tice model with nearest-neighbor interactions: (1) in the 
short-range regime, the critical behavior for the dilute 
model is as without the dilution: it is D = d; (2) in the 
non-mean-field regime it is D G (d : D u ) the system is not 
in the mean-field universality class; (3) in the mean-field 
regime, D > _D U , the system is mean- field. 

Remarkably, this comparison suggests the equality be- 
tween D and the spectral dimension of the Levy graph, 
d. In Ref. 54 it is proved that a fully connected lattice 
with interaction strength decaying as r~~ p (r being the 
distance in a d-dimensional lattice) has spectral dimen- 
sion: 



d = d 

2d 
p — d 



if p> 2 + d 

if p e (d : 2 + d] (6) 



This would prove our conjecture about the three dif- 
ferent critical regimes in the XY model provided that: 
(a) Eq. ^ holds for a Levy lattice with power p, and 
that (b) the relationship between critical properties of a 
model on a graph with spectral dimension d and on a 
lattice of Euclidean dimension d (proved for the Spheri- 
cal and O(oo) modelsP) still holds for the 0(2) model. 
Were these assumptions satisfied, a comparison with Eq. 
H would yield 



(1) 



D = d. 



(7) 



We will later observe, in Sec. |IV| that numerical re- 
sults strongly hints that both assumptions are satisfied. 
In the present section we, rather, provide an analysis 
of the spectral dimension supporting directly the above 
identity, Eq. [7] We numerically estimated the spectral 
dimension of 2D dilute Levy graphs, with several values 
of the power p, through the calculation of the probability 
of self-return of a random walker in the graph after a time 
r, -P(t), a quantity related^ to the spectral dimension 
via 



P(r) 



r -d/2 



(8) 



for large r. Our results are summarized in Fig. [3j in 
which we compare D(p) in Eq. with the estimation 
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FIG. 3: Equivalent short-range dimension -D(p), Eq. |T]), and 
numerically estimated spectral dimension d versus p. The 
light full line is the square lattice case. 



where m is the observable magnetization: 



1 



N 



3=1 



(12) 



and where is a given configuration. Finally, and 

only in the square lattice case, we also measure the second 
moment correlation length ^JSSl 
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(13) 



where k min is the smallest momenta in the Fourier space 
corresponding to the finite-size lattice and x is the 
Fourier transform of the spatial equilibrium two-point 
correlation function 



(14) 



of d at the corresponding p via the histogram of the ran- 
dom walker self-return times. Details of the method are 
reported in App. |Aj 

Fig. [3] suggests the equivalence between the spectral 
dimension of the graph d, and the short-range dimension 
D. This implies that the p-dependence of the spectral 
dimension of the Levy graph with power p is as its equiv- 
alent fully connected version, Eq. [6| and that the 0{2) 
model criticality in a graph is exclusively determined by 
its spectral dimension. 



III. NUMERICAL METHOD, ALGORITHM 
AND DETAILS OF THE SIMULATION 

We now expose the numerical method used to ana- 
lyze the critical properties of the XY model in dilute 2D 
lattices ([3]), via Monte Carlo sampling in finite size real- 
izations of graphs with N vertices and 2N edges. Given 
T and p, we consider both the ensemble average, at T, 
and the graph (topology) average at p: 

(0>=EE °{ S > ex P IS}/? 1 ] V P ( A ( 9 ) 
{J} {S} 

where H is the Hamiltonian of the model Q and O is an 
observable. The following quantities are measured: the 
specific heat 



c = 



1 d(H) 



1 



((H 2 ) (Hf) ; 



(10) 



N dT NT 2 
the susceptibility and the fourth-order Binder cumulant: 



X = A(|m| 2 >, 



Ua 



(|m| 4 ) 
(Imp)* 



(11) 



where i' is such that ry = + r for each i. 

With these observables we analyze the critical prop- 
erties of the model around the critical temperature T c 
using the scaling relations: 



c(T,N) 
X(T,N) 
Ua(T, N) 
&(T,N) 



N't' 9 xitN 1 ' 9 ) 
U i {tN 1 / i ') 



(15) 



where t = T — T c , a and 7 are the standard critical ex- 
ponents and v is the correlation volume exponent. This 
is suited to study scalin g relations in graphs and fully 
connected systemsf 35 * 57 ^ in which the correlation length 
is no longer well-defined, but for which the correlation 
volume V diverges at the critical point as V ~ t~ v . The 
correlation volume exponent is related to the correlation 
length exponent by v = D u v where D u is the upper crit- 
ical dimension of the corresponding short-range system. 

According to the arguments of Sec. [TTJ, and in order 
to elucidate the nature of the phase transition in each 
regime, we have run six sets of simulations with the val- 
ues in the first column of Tab. [nTJ cf. Fig. [2] We 
have studied finite-size realizations of the system in Levy 
graphs with N = L 2 = 2 2n nodes, n = 4, 5, 6, 7, 8. Each 
run, for a fixed topology, consists of 2 21 Monte Carlo 
Steps (MCS). We measure observables each 32 MCS. 
Time averages are performed on exponentially increas- 
ing windows (between 2 k and 2 k+1 , k = 1..., 19,20). 
Topology averages are performed over a sampling of N g 
simulations with different realizations of the graph topol- 
ogy, with N g decreasing for increasing N, going from 
N g = 160 for N = 2 8 to N g = 6 for N = 2 16 . Equilibra- 
tion checks have been done by comparing time averages 
of observables on increasing exponential time windows 
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and verifying the consistency of the energy and magne- 
tization histograms and the correlation length values for 
the last time windows. As a further equilibration test, we 
have checked the coincidence specific heat measurements 
according to both the second and third equalities in Eq. 
[HjJ cf. Fig. top-left panel. 

A. Details of the algorithm 

We now present a more detailed description of the algo- 
rithm used: a home-made high-performance parallel code 
for the Monte Carlo dynamics of spin models defined on 
general networks. The software is developed for CPUs 
architecture and it has been developed with the CUDA 
programming modelP^ A single-spin flip Metropolis up- 
date has been used, with non-connected spins being up- 
dated in parallel by different GPU cores. 

1. Graph coloring 
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FIG. 4: Coloring a dilute Levy 2D graph with power p — 1 
and N = L 2 = 2 8 " 14 nodes with the SLO algorithm. The 
larger the graph, the more homogeneous is the partition of 
the set of graph nodes in subsets corresponding to different 
colors. 



This procedure requires the coloring of each realization 
of the randomly generated graph before dynamics starts. 
The graph nodes are colored with the same color if they 
are not connected to each other. During the simulation, 
sites with a common color are Metropolis-updated syn- 
chronously, and subsets of the set of vertices correspond- 
ing to different colors are processed sequentially on each 
MCS. This is a generalization of the so called Red-Black 
Gauss Seidel algorithm used in the parallelization of spin 
operations in bipartite graphs, as hyper-cubic lattices. 

We approximately color the graph using a var iant of 
the Smallest-Last-Ordering (SLO) algorithm costing 
0[N]. For the simulated sizes (N < 2 16 ) the number 
of colors (equal to two in the p = oo case) turns out 
to be never larger than 6 the worst case (p = 0). As 
one can see m Fig. |4j with our coloring procedure the 
distribution of non-interacting sets becomes more and 
more homogeneous as N increases, thus automatically 
enhancing the algorithm efficiency.^ 
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FIG. 5: Exchange rate between nearby heat baths in the PT 
algorithm in a square lattice with TV = 2 16 sites. The distance 
between consecutive temperatures is AT = 0.005. 



2. Improved equilibrium dynamics 

In our code, besides the Metropolis algorithm, also the 
Parallel Tempering (PT)P and Over-relaxation (OR) 62 
algorithms are implemented. Both algorithms reduce the 
correlation time of the Monte-Carlo Markov processes 
and improve the equilibration.^ 

PT swap attempts are performed (in CPU) every MCS, 
with replicas at different temperatures being updated in 
parallel, as explained in Ref. [331 Fig- [5] illustrates the 
rate of PT swaps between configurations with adjacent 
temperatures at fixed intervals of AT = 0.005, as a func- 
tion of the temperature for a system with N = 2 16 in the 
2D square lattice, for which the critical temperature is 
known to be Tc = 0.8929... . 



3. Memory management 

We have used a storage of the degrees of freedo m in 
the global device memory of the GPU architecture, 1 63 1 64 1 
each thread accessing the 0(2) angle of its corresponding 
graph node in such a way that sites with a common color 
are consecutive in the array, favoring coalesced memory 
access. Each thread, then, accesses an array in global 
memory, from which it reads the list of sites connected to 
the corresponding site. An independent random number 
generator of the Fibonacci type^l is associated to each 
device thread, so that the random numbers are stored in 
shared memory, the values of the random number gen- 
erator parameters being fixed to rp = 250, sp = 128, 
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with which we obtained a satisfactory quality of random 
numbers, as we concluded from the results of our bench- 
mark comparison with the known results for the critical 
temperature of the XY model in the square lattice (see 
the subsection IV A). We used double floating-point pre- 



cision for storing observables, and single precision for the 
calculation of the trigonometric functions in the evalu- 
ation of the energy and magnetization of each site. In 
the latter case we adopted the special fast_math function 
of the GPU architecture, a faster routine specific of the 
GPU architecture.^ 



algorithm 


precision trig 


. function 


time per spin [ns] 


MET + PT 


single 


fast 


1.88 


MET 


single 


fast 


0.635 


MET 


single 


cosf 


0.865 


OR 


single 


fast 


0.36 


OR 


single 


cosf 


0.54 



TABLE II: Computational time of different algorithms per 
spin. It is the total computation time of a run divided by N, 
by the number of copies Nt at different temperatures in the 
PT algorithm and by the number of MCS's. 



4- Computational speed 

We now present some details about the performance of 
our algorithm, referring to a calculation performed in an 
nVidia GPU GTX480 Fermi card. In Tab. Uthe reader 
may find the computation time per spin involved in the 
Metropolis and OR algorithms in a square lattice with 
N = 2 14 sites and PBC, for different choices of the float- 
ing point precision and of the routine used for the com- 
putation of trigonometric operations. In Fig. [6] a com- 
parison of the computation time for the PBC square and 
Levy lattice with p = 7/3 with different sizes is shown. 
Since in a general graph colored with Q colors our algo- 
rithm is nearly Q/2 times slower than the code in the 
square lattice, we also show 2/Q times the computation 
time in the Levy graph for comparison with the square 
lattice case. The minimum time peaks for Metropolis 
and OR algorithms arc 0.55ns and 0.36ns respectively, a 
mark which is competitive with state-of-the-art highly- 
optimized GPU simulations of spin glasses,^ (a direct 
comparison is not possible since their benchmark refers 
to the 0(3) model) and represents a speedup of more that 
five hundred times with respect to a serial C-code running 
on an Intel i7 CPU with 2.67GHz. In the next section 
some accuracy benchmarks of the code are presented. 
All of the simulations in this work have been performed 
using single precision (4 bytes) for the storing of floating- 
point numbers and the fast_math CUDA functions. We 
ran simulations changing both precision and trigonomet- 
ric routines, and introducing OR sweeps, without finding 
essential accuracy improvements. An upgraded and gen- 
eralized version of this algor ithm, designed for the study 
of random laser modes^"-^ in arbitrary topologies, will 
be extensively reviewed and presented in a forthcoming 
workPl 



IV. NUMERICAL RESULTS 

In this section we present numerical results supporting 
the considerations about the XY model criticality as a 
function of p that we presented in section [Tl] All shown 
results are obtained with PBC, since no large quantita- 
tive difference is found adopting FBC (see Sec. [IT|. 
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FIG. 6: Computational time per spin of the Metropolis al- 
gorithm versus N for the square nearest-neighbor lattice and 
p — 7/3 Levy lattice. The serial-CPU run is shown for com- 
parison: a speedup of several hundreds of times can be ob- 
served. 



A. Kosterlitz-Thouless transition: 2D square 
lattice limit 

As a benchmark test we have analyzed the outcome of 
our algorithm in the square lattice (the p — oo limit), 
where the XY model is known to under go an infinitely 
high order transition, the KT transition) 16 * 17 ! The para- 
magnetic high temperature phase, in which vortices un- 
bound, displays exponentially decaying spatial correla- 
tions. The low temperature spin-wave phase is made of 
coupled pairs of vortices of opposite chirality. It is charac- 
terized by the absence of spontaneous magnetization and 
a power-law asymptotic decay of the spatial correlation 
function. Critical exponents depend on the temperature. 

We can compare our results with the analysis of Ref. 
[7T1 finding excellent agreement for all the analyzed quan- 
tities (x, (H), £2, c). As a further test, we have looked at 
several properties of the KT transition, that we will also 
consider as fingerprints for the classification that we will 
do in the following, cf. Sec. |IVE| We enumerate them 
for clarity: 



1. Absence of magnetization. In particular, we 
find numerical evidence of vanishing (|m|)!^l At tem- 
perature larger than the critical point one observes a 
clear behavior (|m|)(T, N) ~ iV -1 / 2 . At lower temper- 
atures, though the decrease with N is slower, it is still 
a decreasing function of N, a behavior compatible with 
(|m|) = in the thermodynamic limit. As we will see 
in the following this behavior is typical for the whole 
short-range regime, cf. Fig. [12} The slow decrease 
is compatible with logarithmic corrections of the type 
X(N) ~ (c + In AO"* + 0[(c + In N)~ 2 }, characteristic of 
the KT transitionPl 

2. Exponential divergence: Kosterlitz-Thouless law. 
Susceptibility and correlation length behave at criticality 
asP 



the relation 



X = X exp 



"A" 

Vt 



(16) 



with t = T - T c > and X = &,X- Figure |§J 
top left panel, illustrates this behavior interpolating the 
largest simulated size. Remarkably, our estimation of 
T c from the fit with the KT law for \ gives T c — 
0.892(9), in excellent agreement with state-of-the-art ac- 
curate estimations^! of T c , =0.8929(4), and the estimate 
from the KT law for £ 2 is T c = 0.90(2). 

3. Continuous specific heat at the transition, as illus- 
trated in Fig. [7J top left panel. 

4. Scale invariance in the whole range T < T c . 

The scaling functions £2 / L and U4 are scale invariant for 
all T < T c in the large- AT limit P^l Differently from what 
happens, for example, in a second-order phase transition, 
one observes a superposition of the finite-size &/L curves 
corresponding to different, sufficiently large, values of L 
(see inset of Fig. [8j top left panel). 



B. Erdos-Renyi limit 

For completeness and checking we have also studied 
the p = limit, corresponding to a random, or Erdos- 
Renyi graph with a Poisson distribution of the degree of 
connectivity and average degree equal to 4. We report 
numerical evidence for a second-order mean-field phase 
transition. Our results are compatible with the theoret- 
ical values of the critical exponents v = 2, 7 = 1 and 
a = 0, as argued in section [TTJ and are in agreement with 
the numerical estimates of Refs. [5F|36( where the mean- 
field transition on Erdos-Renyi graphs of average degree 
3 and 8, respectively, have been analyzed. 

As in the previous subsection, we now present the 
salient analysis for the Erdos-Renyi case. The critical 
temperature can be estimated from the scaling of the 
value at which the finite-size cumulant Ua(T : N) inter- 
sects with C/ 4 (T, 7V/4), a value that we call T(N). As- 
suming the finite-size scaling T(N) = CN- 1 ^ + T c , we 
obtain T c = 1.93(4), in agreement with the analytic value 
T c = 1.9361.^31 The exponent v may be calculated from 



dUi(T,N) 



dT 



oc 



(17) 



C/4— const 



where the derivative of E/4 with respect to the temper- 
ature is evaluated at fixed values of U± in the scaling 
critical region. Approximating the derivative with a fi- 
nite difference (and using the standard deviation of U4 
as an input error for the fit) we obtain v = 2.00(2), in 
agreement with the mean-field value v = 2. Further nu- 
merical evidence for the mean-field values may be found 
in the scaling of the functions C/4, x an d c ( see Fig s -> 
e.g, [9j [Sj [7]). We remark that, although the value of the 
critical exponents is the mean field one, the value of the 
critical temperature is not a universal quantity and does 
not coincide with the Gaussian mean-field value^T = 2. 
We also checked the relation 2/3 + 7 = 2 — a by observing 
the right scaling of the magnetization with the mean- field 
exponent /3 = 1/2. Indeed, an important feature of this 
mean-field transition is that the low-temperature phase 
presents a finite magnetization, and this is confirmed by 
finite size scaling analysis of (|m|) in the low T phase 
(the numerical behaviors are analog to other mean-field 
cases, see, e.g., Fig. 10). 



C. Long-range mean-field regime 

We repeated the same analysis for p = 14/6 and 17/6, 
in the mean-field regime. The first value is nearer to the 
limit of convergence, p — d = 2, of the fully connected 
model, where the largest differences with the dilute model 
could possibly arise. The second value of p is very near 
the mean-field threshold p m f = id/2 = 3. Through finite 
size scaling of the Binder cumulant we estimate T c 
1 .71: 1 and T c = 2.00(1), respectively. The derivative of 
U4 with respect to T allows to estimate v = 2.01(1) and 
v = 1.99(3). These are both compatible with mean-field 
theory. 

Indeed, the system in the mean-field regime behaves 
as a short-range system in D u — 4 dimensions, hence 
displaying the correlation length exponent v = 1/2. The 
correlation volume in the short-range system is the D u -th 
power of the correlation length, and diverges as 



t 



when t — > 0. Since both systems have commo n cr itical 
exponents,^ the above scaling implies, cf. Sec. Ill 



v = — = 2 
v 



(18) 



Let us now consider the correlation length exponent in 
the short-range equivalent system in D(p)-dimensions, 
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FIG. 7: Specific heat versus temperature and t = t iV 1 /" for different Levy lattices and sizes N = 2 2 ™, n = 4, . . . , 8. Top-left: 2D 
PBC XY model. Top-right: dilute model with p = 7/3 in the mean-field regime, v = 2.01(1). Bottom-right: dilute model with 
p — 10/3 in the long-range non-mean-field regime, v = 2.21(1). Bottom-left: dilute model with p = 14/3 in the short-range 
regime. The large N limit of this function is regular and continuous (see the main text). 



r 



that we called v p in section [T] It may be inferred, as in 
reference 1431 comparing the scaling of the free energy den- 
sity of the original system in d dimensions, proportional 
to ST d ~ t~ d / u , and the free energy density of the equiv- 
alent short-range Z?-dimensional system, £~ D ~ t~ D l Vp . 
It follows: 



D(P) 



1 

p — d' 



which coincides with v p found in Ref. 1441 as mentioned 
in Sec. [fl 

In Fig. [7j top-right panel, one observes a scaling of 
the type c(T, N) = c(t ZV 1 / 2 ), suggesting that a = 0. We 
estimated the 7 exponent by fitting it from the numerical 
estimates for the function lnx(T, N) = j / is In N + \n x{t) 
as a function of In TV for a fixed value of t = tN 1 ^ 
in the scaling regime, obtaining 7 = 0.98(4) (c.f. fig- 
ure [8j top-right panel). Finally, one finds (see Fig. 



10 1 that there is spontaneous magnetization in the low- 
temperature phase. 

These numerical results strongly support that the sys- 
tem belongs to the mean-field universality class in the 
range p € [0 : 3<i/2]. In the whole range the critical 
exponents 9 = 2, 7=1, a = are well defined. On 
the other hand, the correlation length exponent of the 
equivalent short-range system, v p in Eq. 



19 is no longer 



(19) defined for p < d, where the dimensionality D diverges 



D. Long-range non- mean-field regime 

An analogous analysis leads to a different behavior in 
the non-mean-field regime, for p = 10/3. For the critical 
temperature we estimate T c = 1.749(6). The finite-size 
scaling analysis of Eq. (17) reveals a correlation volume 
exponent larger that the mean-field value 2: D = 2.21(1) 
(see Fig. 11). Estimating the 7 exponent from the N 
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FIG. 8: Magnetic susceptibility versus temperature for different Levy graphs (p and L values are as in Fig. [7}. Top-left: 2D 
PBC XY model. As the size increases, x{T) becomes more and more similar to the functional law Eq. For N = 2 16 it 

is lnxo = —1-48(2), b x — 2.89(1) and T c — 0.8929(4). Inset: correlation length. Also in this case the large size limit obeys 
Eq. (16|, with ln£o = —6.5(7) and b% = 1.4(4). In the low-temperature phase the £2 curves collapse. Top-right: Scaling of 
the susceptibility for the mean-field Levy model at p — 7/3, cf . Eq. (15 1. 

to 7 = 0.98(4). Bottom-right: Scaling of the susceptibility for the long-range non-mean-field model at p 



Inset: finite size behavior of lnx at fixed t, leading 

10/3 for T > T c . 



The critical exponents are v = 2.21(1), 7 = 1.34(2). Bottom-left: susceptibility vs. temperature for the Levy short-range case, 
p = 14/3. The KT fit gives lnxo = -1.52(2), b x = 2.80(2) and T c = 1.338(1). 



dependence of the function In x[T, N) for a fixed value of 
the scaling variable t, we obtain 7 = 1.34(1), cf., Fig. § 
bottom-right panel, different from the mean-field value, 
1. On the other hand, the specific heat data is compat- 
ible with a — 0. Indeed, although the TV-dependence 
of the specific heat is slower with respect to the mean- 
field case, the data in Fig. [7J bottom-right panel, suggest 
that curves corresponding to different iV's collapse in the 
large- A^ limit for common values of t. Finally, we observe 
a low-T magnetized state, (|m|) 7^ 0. 



E. Short-range Regime 

The case p = 14/3 in the short-range regime presents 
all the features that we presented as numerical finger- 



prints of the KT transition in subsection IV A (1) the 
low-T phase is unmagnetized, as deduced from the de- 
creasing behavior of (|m|) vs. A^ for the sizes studied, cf. 
Fig. [l2j (2) The KT law, Eq. ^ is satisfied, yielding 
an estimate T c = 1.338(1), cf. Fig. [8j bottom-left panel; 
(3) the specific heat is not divergent nor discontinuous at 
the transition, cf., Fig. [7J bottom- left panel, at difference 
with the other cases, in which the discontinuity of c was 
reflected in the discontinuity of the derivative of its scal- 
ing function, c; (4) A- numerical evidence of the fact that 
the quantity U4 is scale invariant for sufficiently large 
sizes is shown in Fig. [T3j in which one observes that the 
quantity \U^{T,N) — U±(T, AN)\ decreases for increasing 
values of A^ at constant values of T immediately below 
the transition. 

We summarize all our results in Tab. IIIII 
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T 

FIG. 12: Magnetization versus temperature for p = 14/3 in 
FIG. 10: Magnetization versus temperature for p = 7/3 in the short-range regime. In the low-T phase the magnetization 
the mean-field regime. The low temperature phase exhibits monotonously decreases with N, indicating that this phase 
spontaneous magnetization. does not exhibit spontaneous magnetization. 



V. CONCLUSIONS AND PERSPECTIVES 

Extensive numerical simulations on the 2D Levy lat- 
tice yield evidence for three different critical regimes, 
corresponding to given intervals in the Levy exponent 
p governing the topology of the graphs: short-range 
p £ [2+d, oo), non-mean-field long-range p € (3d/2, 2+d) 
and mean-field p e [0, 3d/ 2]. For each value of p, the crit- 
ical behavior is in a one-to-one correspondence with that 
of a short-range XY model in D dimensions. This short- 
range effective dimension is D = d for p > 2 + d, while it 
is D = 2d/(p — d) for p e (d, 2 + d), in the long-range and 
part of the mean-field regime. As p — > d, it tends to infi- 
nite. This is the value of p for which the fully connected 



version of the model displays a divergent energy. 

In the mean-field regime we measured the standard 
mean-field exponents 7 = 1, a = and 9 = 2. The 
latter is the correlation volume exponent, related with 
the correlation length exponent v = v/D n = 1/2. These 
exponents agree with the already mentioned theoretical 
predictions for the D-dimensional equivalent model in 



the mean-field regime: v p , rj p , 7 P (see subsection IV C I. 
In the non-mean field regime, instead, we find a continu- 
ous phase transition with different critical exponents and 
a low-temperature phase exhibiting spontaneous symme- 
try breaking. Finally, we report evidences for the onset 
of a KT-like transition in the short-range regime. 

The hypothesis presented in Sec. In) that D is the 
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understanding of the physics of interacting systems on 
non-regular structures and extending the known univer- 
sal properties of the spectral dimension of the 0(n — > oo) 
model^ to finite n. 

Summarizing, we have found that the critical proper- 
ties of the 0(2) model in a graph can be divided in three 
regimes characterized by the spectral dimension of the 
graph, which plays the role of the dimension of a short- 
range model with common critical properties. As a per- 
spective, we propose to investigate how the introduction 
of disorder changes this scenario.^ This is the object of 
an ongoing research. 



T 
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FIG. 13: Binder cumulant versus temperature for p — 14/3 
in the short-range regime. There is evidence for the scale- 
invariance of this quantity for T < T c in the large- N limit. 



p 





14/6 


17/6 


20/6 


28/6 


oo 


T c 


1.932(1) 


1.94(1) 


2.00(1) 


1.749(6) 


1.338(1) 


0.892(9) 


V 


2.00(2) 


2.01(1) 


1.99(3) 


2.21(1) 






7 


1.00(1) 


0.98(4) 


1.04(8) 


1.34(2) 






a 











< 0.1 






(|m|> 




t^O 






= 


= 


UC 


MF 


MF 


MF 


NMF 


SR 


SR 



TABLE III: Numerical results for the critical exponents and 
the critical properties of the XY model in different Levy 
graphs, p = oo stays for the 2D square lattice with PBC. UC 
stays for "universality class": MF=mean field, NMF=non- 
mean-field, SR=short-range. 



spectral dimension of the dilute graph, appears to be 
confirmed in the three regimes. The numerical data pre- 



sented in section IV suggest that the low-temperature 
phase of the XY Levy model in the short-range regime 
is not magnetized, while there is a spontaneous symme- 
try breaking for p < p sr = 2 + d. The already mentioned 
gener alized Mermin- Wagner theorem for the XY model in 
graph s 14 * 15 ^ would, then, imply that d = 2 for p > 4, while 
d > 2 for p < 4 in the dilute lattice, in agreement with 
Eq. M. On the other hand, for p < 4, the identification 
D = a appears to be confirmed by a numerical estima- 
tion, see figure[3]and App. [A] Furthermore, for p — » d, D 
diverges. An infinite spectral dimension, indeed, occurs 
in the Bethe lattice limit and, generally, in any graphs 
not satisfying the polynomial growth condition.ES 

The result D = d would imply that the critical be- 
havior of the 0(2) model on a graph characterized by a 
spectral dimension d coincides with the one of the short- 
range d-dimensional 0(2) model, allowing for a deeper 
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Appendix A: Spectral density estimation 

In order to estimate d we have studied the histogram 
of return times of a collectivity of 4.10 4 random walkers 
averaging over 20 different finite size realizations of 2D 
dilute graphs with sizes N = 128 2 and 196 2 and 6 dif- 
ferent values of the the power p = oo, 5, 4.5, 4, 3.8 and 
11/3. Such a histogram is proportional to the probabil- 
ity of self-return of a random walker in the graph after 
a time r: P(r). For large values of p > 2 + d and in 
the square lattice case the estimation of d is very accu- 
rate since the function P presents a very clear power-law 
behavior even at large times and finite-size effects are 
not an issue. For smaller values of p, however, such a 
measure becomes less and less accurate. This is due to 
the presence of "shortcuts" on the graph that take the 
walker to the boundaries of the original lattice, where 
the probability P is overestimated, and to the existence 
of low-connected nodes. To cope with these drawbacks, 
for p > 4 our random walkers start from the center of the 
original finite-size 2D lattice, while for p < 4, each real- 
ization of the walker starts from a node which is chosen 
at random between the subset of nodes whose degree is 
larger than two. In Fig. [14] we present the P histograms 
(up to an arbitrary, p-dependent constant) for each stud- 
ied value of p, from which the data of Fig. [3] have been 
inferred. 
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FIG. 14: Logarithm of the histogram of self-return times of 
a random walker in dilute 2D graphs of size N = 196 2 with 
different values of p. Points are the numerical measures, while 
lines are the fits with a linear function (their endpoints indi- 
cate the corresponding fit intervals). The slope of the fits is 
— 2 times the spectral dimension (see figure [3b. 
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